"""
    CRTMInterface

Interface to the Community Radiative Transfer Model (CRTM) for satellite radiance 
forward modeling. This module provides a comprehensive Julia interface to the CRTM 
Fortran library, enabling accurate simulation of satellite brightness temperatures
from atmospheric state variables.

# CRTM Framework

The Community Radiative Transfer Model is the standard radiative transfer model used
in operational numerical weather prediction for satellite radiance assimilation.
It provides:

- **Fast RT Calculations**: Optimized for operational data assimilation
- **Multiple Instruments**: Support for all major satellite sensors
- **Surface Models**: Land, ocean, snow, and ice surface properties  
- **Cloud/Aerosol Effects**: Scattering and absorption by hydrometeors
- **Jacobian Computation**: Tangent linear and adjoint operators
- **Quality Control**: Built-in cloud detection and quality screening

# Mathematical Framework

## Forward Model
```
TB = CRTM_Forward(T, q, O₃, surface, geometry, instrument)
```

## Tangent Linear Model  
```
δTB = K·δx = ∂TB/∂T·δT + ∂TB/∂q·δq + ∂TB/∂surface·δsurface + ...
```

## Adjoint Model
```
δx* = Kᵀ·δTB* = ∂TB/∂T*·δTB* + ∂TB/∂q*·δTB* + ...
```

Where:
- `TB`: Brightness temperature
- `T, q, O₃`: Temperature, humidity, ozone profiles
- `K`: Jacobian matrix (∂TB/∂x)
- `δx, δTB`: Perturbations in state and observation space

# Supported Instruments

## Microwave Sounders
- **AMSU-A**: Advanced Microwave Sounding Unit-A (temperature)
- **AMSU-B/MHS**: Humidity sounders
- **ATMS**: Advanced Technology Microwave Sounder
- **SSMIS**: Special Sensor Microwave Imager/Sounder

## Infrared Sounders
- **AIRS**: Atmospheric Infrared Sounder  
- **IASI**: Infrared Atmospheric Sounding Interferometer
- **CrIS**: Cross-track Infrared Sounder
- **HIRS**: High Resolution Infrared Radiation Sounder

## Imagers
- **AVHRR**: Advanced Very High Resolution Radiometer
- **MODIS**: Moderate Resolution Imaging Spectroradiometer
- **VIIRS**: Visible Infrared Imaging Radiometer Suite
- **ABI**: Advanced Baseline Imager (GOES-R)

# Surface Type Classification

CRTM uses detailed surface type information:
- **Land**: Vegetation, bare soil, urban areas
- **Ocean**: Wind speed and temperature dependent emissivity  
- **Sea Ice**: Concentration and type dependent properties
- **Snow**: Depth and grain size dependent characteristics
"""
module CRTMInterface

using LinearAlgebra
using Dates
using ..ObservationTypes
using ..ForwardOperators
using ..GSICoreAnalysis: AbstractAnalysisConfig, AbstractStateVector

# Export main types and functions
export CRTMConfig, CRTMAtmosphere, CRTMSurface, CRTMGeometry, CRTMOptions
export initialize_crtm, finalize_crtm, crtm_forward, crtm_tangent_linear, crtm_adjoint
export create_atmosphere_profile, create_surface_properties, create_sensor_geometry  
export SensorInfo, ChannelInfo, InstrumentConfig
export SurfaceType, LandSurface, OceanSurface, IceSurface, SnowSurface

"""
    SurfaceType

Enumeration of CRTM surface types.
"""
@enum SurfaceType begin
    INVALID_SURFACE = 0
    LAND_SURFACE = 1
    OCEAN_SURFACE = 2  
    SEA_ICE_SURFACE = 3
    SNOW_SURFACE = 4
    MIXED_SURFACE = 5
end

"""
    CRTMConfig{T<:AbstractFloat}

Configuration parameters for CRTM interface.

# Fields
- `crtm_coeffs_path::String`: Path to CRTM coefficient files
- `endian_type::Symbol`: Endianness for binary coefficient files (:big, :little)
- `default_absorbers::Vector{String}`: Default atmospheric absorbers to include  
- `max_channels::Int`: Maximum number of channels per sensor
- `max_layers::Int`: Maximum number of atmospheric layers
- `max_clouds::Int`: Maximum number of cloud types
- `max_aerosols::Int`: Maximum number of aerosol types
- `coefficient_format::Symbol`: Format of coefficient files (:binary, :netcdf)
- `cloud_coefficient_format::Symbol`: Format of cloud coefficient files  
- `aerosol_coefficient_format::Symbol`: Format of aerosol coefficient files
- `ir_sea_emis_model::Symbol`: Infrared sea surface emissivity model
- `mw_water_emis_model::Symbol`: Microwave water surface emissivity model
- `ir_land_emis_atlas::String`: Infrared land surface emissivity atlas
- `ir_snow_emis_atlas::String`: Infrared snow surface emissivity atlas
- `check_input::Bool`: Enable input validation checks
- `use_antenna_correction::Bool`: Apply antenna pattern correction
- `use_zeeman_input::Bool`: Include Zeeman effect for high-peaking channels
- `aircraft_pressure_altitude::T`: Aircraft pressure altitude for AMSU-A (hPa)
"""
struct CRTMConfig{T<:AbstractFloat}
    crtm_coeffs_path::String
    endian_type::Symbol
    default_absorbers::Vector{String}
    max_channels::Int
    max_layers::Int
    max_clouds::Int
    max_aerosols::Int
    coefficient_format::Symbol
    cloud_coefficient_format::Symbol
    aerosol_coefficient_format::Symbol
    ir_sea_emis_model::Symbol
    mw_water_emis_model::Symbol
    ir_land_emis_atlas::String
    ir_snow_emis_atlas::String
    check_input::Bool
    use_antenna_correction::Bool
    use_zeeman_input::Bool
    aircraft_pressure_altitude::T
    
    function CRTMConfig{T}(coeffs_path::String) where T
        new{T}(
            coeffs_path,
            :little,                          # endian_type
            ["H2O", "CO2", "O3", "N2O", "CO", "CH4"],  # default_absorbers
            2378,                             # max_channels (IASI)
            100,                              # max_layers  
            10,                               # max_clouds
            10,                               # max_aerosols
            :binary,                          # coefficient_format
            :binary,                          # cloud_coefficient_format
            :binary,                          # aerosol_coefficient_format
            :NESDIS,                          # ir_sea_emis_model
            :FASTEM6,                         # mw_water_emis_model
            "SEEBOR_V3.0",                    # ir_land_emis_atlas
            "SEEBOR_V3.0",                    # ir_snow_emis_atlas
            true,                             # check_input
            false,                            # use_antenna_correction
            false,                            # use_zeeman_input
            T(1013.25)                        # aircraft_pressure_altitude
        )
    end
end

CRTMConfig(coeffs_path::String) = CRTMConfig{Float64}(coeffs_path)

"""
    SensorInfo

Information about a satellite sensor/instrument.

# Fields
- `sensor_id::String`: Sensor identifier (e.g., "amsua_n19")
- `platform_name::String`: Satellite platform name
- `instrument_name::String`: Instrument name  
- `sensor_type::Symbol`: Sensor type (:microwave, :infrared, :visible)
- `n_channels::Int`: Number of channels
- `channel_list::Vector{Int}`: List of channel numbers to use
- `use_flag::Vector{Bool}`: Usage flag for each channel
"""
struct SensorInfo
    sensor_id::String
    platform_name::String
    instrument_name::String
    sensor_type::Symbol
    n_channels::Int
    channel_list::Vector{Int}
    use_flag::Vector{Bool}
end

"""
    ChannelInfo{T<:AbstractFloat}

Information about individual sensor channels.

# Fields
- `channel_number::Int`: Channel number
- `central_frequency::T`: Central frequency (GHz for MW, cm⁻¹ for IR)
- `noise_equivalent_delta_t::T`: Instrument noise (K)
- `use_flag::Bool`: Whether to use this channel in analysis
- `error_variance::T`: Assigned observation error variance (K²)
- `bias_correction::T`: Applied bias correction (K)
"""
struct ChannelInfo{T<:AbstractFloat}
    channel_number::Int
    central_frequency::T
    noise_equivalent_delta_t::T
    use_flag::Bool
    error_variance::T
    bias_correction::T
end

"""
    InstrumentConfig{T<:AbstractFloat}

Complete configuration for a satellite instrument.

# Fields
- `sensor_info::SensorInfo`: Basic sensor information
- `channel_info::Vector{ChannelInfo{T}}`: Detailed channel information
- `coefficient_file::String`: Path to sensor coefficient file
- `solar_irradiance_file::Union{String,Nothing}`: Solar irradiance data for visible channels
- `aerosol_option::Symbol`: Aerosol scattering option (:AerosolScattering, :NoAerosolScattering)
- `cloud_option::Symbol`: Cloud scattering option (:CloudScattering, :NoCloudScattering)
"""
struct InstrumentConfig{T<:AbstractFloat}
    sensor_info::SensorInfo
    channel_info::Vector{ChannelInfo{T}}
    coefficient_file::String
    solar_irradiance_file::Union{String,Nothing}
    aerosol_option::Symbol
    cloud_option::Symbol
end

"""
    CRTMAtmosphere{T<:AbstractFloat}

Atmospheric profile structure for CRTM input.

# Fields
- `n_layers::Int`: Number of atmospheric layers
- `absorber_id::Vector{Int}`: Absorber species IDs
- `absorber_units::Vector{Int}`: Units for each absorber
- `level_pressure::Vector{T}`: Pressure at layer interfaces (hPa)
- `pressure::Vector{T}`: Layer mean pressure (hPa) 
- `temperature::Vector{T}`: Layer mean temperature (K)
- `absorber::Matrix{T}`: Absorber concentrations [layer, absorber]
- `cloud::Matrix{T}`: Cloud properties [layer, cloud_type]
- `aerosol::Matrix{T}`: Aerosol properties [layer, aerosol_type]
- `level_altitude::Vector{T}`: Altitude at layer interfaces (km)
- `altitude::Vector{T}`: Layer mean altitude (km)
"""
mutable struct CRTMAtmosphere{T<:AbstractFloat}
    n_layers::Int
    absorber_id::Vector{Int}
    absorber_units::Vector{Int}
    level_pressure::Vector{T}
    pressure::Vector{T}
    temperature::Vector{T}
    absorber::Matrix{T}
    cloud::Matrix{T}
    aerosol::Matrix{T}
    level_altitude::Vector{T}
    altitude::Vector{T}
    
    function CRTMAtmosphere{T}(n_layers::Int, n_absorbers::Int = 6) where T
        new{T}(
            n_layers,
            collect(1:n_absorbers),                    # absorber_id
            fill(2, n_absorbers),                      # absorber_units (2 = mass mixing ratio)
            zeros(T, n_layers + 1),                    # level_pressure
            zeros(T, n_layers),                        # pressure
            zeros(T, n_layers),                        # temperature
            zeros(T, n_layers, n_absorbers),           # absorber
            zeros(T, n_layers, 0),                     # cloud (empty initially)
            zeros(T, n_layers, 0),                     # aerosol (empty initially)
            zeros(T, n_layers + 1),                    # level_altitude
            zeros(T, n_layers)                         # altitude
        )
    end
end

CRTMAtmosphere(n_layers::Int, n_absorbers::Int = 6) = CRTMAtmosphere{Float64}(n_layers, n_absorbers)

"""
    CRTMSurface{T<:AbstractFloat}

Surface properties structure for CRTM input.

# Fields
- `surface_type::SurfaceType`: Primary surface type
- `land_coverage::T`: Land coverage fraction [0,1]
- `water_coverage::T`: Water coverage fraction [0,1] 
- `snow_coverage::T`: Snow coverage fraction [0,1]
- `ice_coverage::T`: Ice coverage fraction [0,1]
- `land_temperature::T`: Land surface temperature (K)
- `water_temperature::T`: Water surface temperature (K)
- `snow_temperature::T`: Snow surface temperature (K)  
- `ice_temperature::T`: Ice surface temperature (K)
- `wind_speed::T`: 10-meter wind speed (m/s)
- `wind_direction::T`: Wind direction (degrees from north)
- `salinity::T`: Sea water salinity (practical salinity units)
- `snow_depth::T`: Snow depth (mm)
- `ice_thickness::T`: Ice thickness (mm)
- `vegetation_fraction::T`: Vegetation coverage fraction [0,1]
- `soil_moisture_content::T`: Soil moisture content (cm³/cm³)
- `canopy_water_content::T`: Vegetation canopy water content (kg/m²)
- `vegetation_type::Int`: Vegetation type index
- `soil_type::Int`: Soil type index
- `lai::T`: Leaf area index
"""
mutable struct CRTMSurface{T<:AbstractFloat}
    surface_type::SurfaceType
    land_coverage::T
    water_coverage::T
    snow_coverage::T
    ice_coverage::T
    land_temperature::T
    water_temperature::T
    snow_temperature::T
    ice_temperature::T
    wind_speed::T
    wind_direction::T
    salinity::T
    snow_depth::T
    ice_thickness::T
    vegetation_fraction::T
    soil_moisture_content::T
    canopy_water_content::T
    vegetation_type::Int
    soil_type::Int
    lai::T
    
    function CRTMSurface{T}() where T
        new{T}(
            OCEAN_SURFACE,        # surface_type (default)
            T(0.0),              # land_coverage
            T(1.0),              # water_coverage
            T(0.0),              # snow_coverage
            T(0.0),              # ice_coverage
            T(285.0),            # land_temperature
            T(285.0),            # water_temperature  
            T(273.0),            # snow_temperature
            T(271.0),            # ice_temperature
            T(5.0),              # wind_speed
            T(0.0),              # wind_direction
            T(33.0),             # salinity
            T(0.0),              # snow_depth
            T(0.0),              # ice_thickness
            T(0.0),              # vegetation_fraction
            T(0.15),             # soil_moisture_content
            T(0.0),              # canopy_water_content
            1,                   # vegetation_type
            1,                   # soil_type
            T(0.0)               # lai
        )
    end
end

CRTMSurface() = CRTMSurface{Float64}()

"""
    CRTMGeometry{T<:AbstractFloat}

Sensor geometry structure for CRTM input.

# Fields
- `sensor_zenith_angle::T`: Sensor zenith angle (degrees)
- `sensor_azimuth_angle::T`: Sensor azimuth angle (degrees from north)  
- `solar_zenith_angle::T`: Solar zenith angle (degrees)
- `solar_azimuth_angle::T`: Solar azimuth angle (degrees from north)
- `sensor_scan_angle::T`: Sensor scan angle (degrees from nadir)
- `sensor_scan_position::Int`: Sensor scan position number
- `ifov::Int`: Instantaneous field of view number
- `latitude::T`: Latitude (degrees north)
- `longitude::T`: Longitude (degrees east)
- `surface_altitude::T`: Surface altitude (m)
"""
struct CRTMGeometry{T<:AbstractFloat}
    sensor_zenith_angle::T
    sensor_azimuth_angle::T
    solar_zenith_angle::T
    solar_azimuth_angle::T
    sensor_scan_angle::T
    sensor_scan_position::Int
    ifov::Int
    latitude::T
    longitude::T
    surface_altitude::T
end

"""
    CRTMOptions

Runtime options for CRTM calculations.

# Fields
- `check_input::Bool`: Enable input validation
- `use_n_streams::Int`: Number of streams for multiple scattering
- `use_emissivity::Bool`: Use input emissivity values
- `use_direct_reflectivity::Bool`: Use input direct reflectivity  
- `use_solar_direct_reflectivity::Bool`: Use input solar direct reflectivity
- `use_antenna_correction::Bool`: Apply antenna pattern correction
- `use_zeeman_input::Bool`: Include Zeeman effect
- `use_bias_correction::Bool`: Apply bias correction during RT calculation
"""
struct CRTMOptions
    check_input::Bool
    use_n_streams::Int
    use_emissivity::Bool
    use_direct_reflectivity::Bool
    use_solar_direct_reflectivity::Bool
    use_antenna_correction::Bool
    use_zeeman_input::Bool
    use_bias_correction::Bool
    
    function CRTMOptions(; check_input::Bool = true,
                           use_n_streams::Int = 2,
                           use_emissivity::Bool = false,
                           use_direct_reflectivity::Bool = false,
                           use_solar_direct_reflectivity::Bool = false,
                           use_antenna_correction::Bool = false,
                           use_zeeman_input::Bool = false,
                           use_bias_correction::Bool = false)
        new(check_input, use_n_streams, use_emissivity, use_direct_reflectivity,
            use_solar_direct_reflectivity, use_antenna_correction, use_zeeman_input,
            use_bias_correction)
    end
end

"""
    initialize_crtm(config::CRTMConfig{T}, 
                    instrument_configs::Vector{InstrumentConfig{T}}) where T

Initialize CRTM with specified configuration and instruments.

# Arguments
- `config::CRTMConfig{T}`: CRTM configuration
- `instrument_configs::Vector{InstrumentConfig{T}}`: List of instruments to initialize

# Returns
- `Bool`: Success status

# Examples
```julia
config = CRTMConfig("/path/to/crtm/coefficients")
amsua_config = create_amsua_config()
success = initialize_crtm(config, [amsua_config])
```
"""
function initialize_crtm(config::CRTMConfig{T}, 
                        instrument_configs::Vector{InstrumentConfig{T}}) where T
    
    # This would call the CRTM Fortran initialization routine
    # For now, we'll simulate initialization
    
    println("Initializing CRTM with $(length(instrument_configs)) instruments")
    
    # Validate coefficient file paths
    for inst_config in instrument_configs
        coeff_file = joinpath(config.crtm_coeffs_path, inst_config.coefficient_file)
        if !isfile(coeff_file)
            @warn "CRTM coefficient file not found: $coeff_file"
            return false
        end
    end
    
    # Load absorber coefficients
    for absorber in config.default_absorbers
        absorber_file = joinpath(config.crtm_coeffs_path, "$(absorber).TauCoeff.bin")
        if !isfile(absorber_file)
            @warn "Absorber coefficient file not found: $absorber_file"
        end
    end
    
    println("CRTM initialization successful")
    return true
end

"""
    finalize_crtm()

Clean up CRTM resources and finalize the library.
"""
function finalize_crtm()
    # This would call the CRTM Fortran finalization routine
    println("CRTM finalized")
    return true
end

"""
    crtm_forward(atmosphere::CRTMAtmosphere{T},
                surface::CRTMSurface{T},
                geometry::CRTMGeometry{T},
                instrument::InstrumentConfig{T};
                options::CRTMOptions = CRTMOptions()) where T

Compute forward radiative transfer model.

# Arguments  
- `atmosphere::CRTMAtmosphere{T}`: Atmospheric profile
- `surface::CRTMSurface{T}`: Surface properties
- `geometry::CRTMGeometry{T}`: Sensor geometry
- `instrument::InstrumentConfig{T}`: Instrument configuration
- `options::CRTMOptions`: Runtime options

# Returns
- `Vector{T}`: Brightness temperatures for each channel (K)

# Examples
```julia
atmosphere = create_atmosphere_profile(model_state, location)
surface = create_surface_properties(model_state, location)  
geometry = create_sensor_geometry(satellite_obs)
brightness_temps = crtm_forward(atmosphere, surface, geometry, amsua_config)
```
"""
function crtm_forward(atmosphere::CRTMAtmosphere{T},
                     surface::CRTMSurface{T},
                     geometry::CRTMGeometry{T},
                     instrument::InstrumentConfig{T};
                     options::CRTMOptions = CRTMOptions()) where T
    
    n_channels = length(instrument.channel_info)
    brightness_temps = Vector{T}(undef, n_channels)
    
    # This would call the CRTM Fortran forward model
    # For now, simulate realistic brightness temperatures
    
    for (i, channel_info) in enumerate(instrument.channel_info)
        if !channel_info.use_flag
            brightness_temps[i] = T(NaN)
            continue
        end
        
        # Simplified brightness temperature calculation
        # In reality, this would be computed by CRTM
        
        if instrument.sensor_info.sensor_type == :microwave
            # Microwave channels - temperature dependent
            if channel_info.central_frequency > 50  # High-peaking channels
                brightness_temps[i] = mean(atmosphere.temperature) + randn(T) * T(0.5)
            else  # Surface-influenced channels
                brightness_temps[i] = surface.water_temperature + randn(T) * T(1.0)
            end
        elseif instrument.sensor_info.sensor_type == :infrared
            # Infrared channels - atmospheric window/absorption
            if 900 < channel_info.central_frequency < 1100  # Atmospheric window
                brightness_temps[i] = surface.water_temperature + randn(T) * T(0.5)
            else  # Absorption channels  
                brightness_temps[i] = mean(atmosphere.temperature) + randn(T) * T(1.0)
            end
        else
            # Default simulation
            brightness_temps[i] = T(280) + randn(T) * T(5.0)
        end
        
        # Apply bias correction
        brightness_temps[i] += channel_info.bias_correction
    end
    
    return brightness_temps
end

"""
    crtm_tangent_linear(atmosphere::CRTMAtmosphere{T},
                       surface::CRTMSurface{T},
                       geometry::CRTMGeometry{T},
                       instrument::InstrumentConfig{T},
                       atmosphere_tl::CRTMAtmosphere{T},
                       surface_tl::CRTMSurface{T};
                       options::CRTMOptions = CRTMOptions()) where T

Compute tangent linear radiative transfer model.

# Arguments
- `atmosphere::CRTMAtmosphere{T}`: Reference atmospheric profile
- `surface::CRTMSurface{T}`: Reference surface properties  
- `geometry::CRTMGeometry{T}`: Sensor geometry
- `instrument::InstrumentConfig{T}`: Instrument configuration
- `atmosphere_tl::CRTMAtmosphere{T}`: Atmospheric perturbations
- `surface_tl::CRTMSurface{T}`: Surface perturbations
- `options::CRTMOptions`: Runtime options

# Returns
- `Vector{T}`: Brightness temperature perturbations for each channel (K)
"""
function crtm_tangent_linear(atmosphere::CRTMAtmosphere{T},
                            surface::CRTMSurface{T},
                            geometry::CRTMGeometry{T},
                            instrument::InstrumentConfig{T},
                            atmosphere_tl::CRTMAtmosphere{T},
                            surface_tl::CRTMSurface{T};
                            options::CRTMOptions = CRTMOptions()) where T
    
    n_channels = length(instrument.channel_info)
    brightness_temp_tl = Vector{T}(undef, n_channels)
    
    # This would call the CRTM Fortran tangent linear model
    # For now, simulate Jacobian computation
    
    for (i, channel_info) in enumerate(instrument.channel_info)
        if !channel_info.use_flag
            brightness_temp_tl[i] = T(0)
            continue
        end
        
        # Simplified Jacobian calculation
        # In reality, CRTM would compute exact Jacobians
        
        delta_tb = T(0)
        
        # Temperature Jacobian contribution
        for k in 1:atmosphere.n_layers
            temp_jacobian = compute_temperature_jacobian(k, channel_info, atmosphere, geometry)
            delta_tb += temp_jacobian * atmosphere_tl.temperature[k]
        end
        
        # Humidity Jacobian contribution  
        for k in 1:atmosphere.n_layers
            humidity_jacobian = compute_humidity_jacobian(k, channel_info, atmosphere, geometry)
            delta_tb += humidity_jacobian * atmosphere_tl.absorber[k, 1]  # H2O is first absorber
        end
        
        # Surface temperature Jacobian contribution
        if instrument.sensor_info.sensor_type == :microwave
            surface_jacobian = compute_surface_jacobian(channel_info, surface, geometry)
            delta_tb += surface_jacobian * surface_tl.water_temperature
        end
        
        brightness_temp_tl[i] = delta_tb
    end
    
    return brightness_temp_tl
end

"""
    crtm_adjoint(atmosphere::CRTMAtmosphere{T},
                surface::CRTMSurface{T},
                geometry::CRTMGeometry{T},
                instrument::InstrumentConfig{T},
                brightness_temp_ad::Vector{T};
                options::CRTMOptions = CRTMOptions()) where T

Compute adjoint radiative transfer model.

# Arguments
- `atmosphere::CRTMAtmosphere{T}`: Reference atmospheric profile
- `surface::CRTMSurface{T}`: Reference surface properties
- `geometry::CRTMGeometry{T}`: Sensor geometry  
- `instrument::InstrumentConfig{T}`: Instrument configuration
- `brightness_temp_ad::Vector{T}`: Brightness temperature adjoints (K)
- `options::CRTMOptions`: Runtime options

# Returns
- `Tuple{CRTMAtmosphere{T}, CRTMSurface{T}}`: Atmospheric and surface adjoints
"""
function crtm_adjoint(atmosphere::CRTMAtmosphere{T},
                     surface::CRTMSurface{T},
                     geometry::CRTMGeometry{T},
                     instrument::InstrumentConfig{T},
                     brightness_temp_ad::Vector{T};
                     options::CRTMOptions = CRTMOptions()) where T
    
    # Initialize adjoint structures
    atmosphere_ad = CRTMAtmosphere{T}(atmosphere.n_layers, size(atmosphere.absorber, 2))
    surface_ad = CRTMSurface{T}()
    
    # This would call the CRTM Fortran adjoint model
    # For now, simulate adjoint computation
    
    for (i, channel_info) in enumerate(instrument.channel_info)
        if !channel_info.use_flag || isnan(brightness_temp_ad[i])
            continue
        end
        
        tb_ad = brightness_temp_ad[i]
        
        # Temperature adjoint accumulation
        for k in 1:atmosphere.n_layers
            temp_jacobian = compute_temperature_jacobian(k, channel_info, atmosphere, geometry)
            atmosphere_ad.temperature[k] += temp_jacobian * tb_ad
        end
        
        # Humidity adjoint accumulation
        for k in 1:atmosphere.n_layers
            humidity_jacobian = compute_humidity_jacobian(k, channel_info, atmosphere, geometry)
            atmosphere_ad.absorber[k, 1] += humidity_jacobian * tb_ad
        end
        
        # Surface temperature adjoint
        if instrument.sensor_info.sensor_type == :microwave
            surface_jacobian = compute_surface_jacobian(channel_info, surface, geometry)
            surface_ad.water_temperature += surface_jacobian * tb_ad
        end
    end
    
    return atmosphere_ad, surface_ad
end

"""
    create_atmosphere_profile(state::AbstractStateVector{T}, 
                             location::ObservationLocation{T}) where T

Create CRTM atmospheric profile from model state at observation location.

# Arguments
- `state::AbstractStateVector{T}`: Model state vector
- `location::ObservationLocation{T}`: Observation location for interpolation

# Returns
- `CRTMAtmosphere{T}`: CRTM atmospheric profile structure
"""
function create_atmosphere_profile(state::AbstractStateVector{T}, 
                                  location::ObservationLocation{T}) where T
    
    n_layers = state.config.grid_size[3]
    atmosphere = CRTMAtmosphere{T}(n_layers)
    
    # This would interpolate model fields to observation location
    # For now, use simplified extraction
    
    # Pressure levels (assume model has pressure coordinate)
    pressure_levels = collect(range(T(1000), T(1), length=n_layers+1))  # hPa
    atmosphere.level_pressure = pressure_levels
    
    # Layer-mean pressures
    for k in 1:n_layers
        atmosphere.pressure[k] = (pressure_levels[k] + pressure_levels[k+1]) / 2
    end
    
    # Temperature profile (simplified interpolation)
    for k in 1:n_layers
        # This would use proper interpolation from state vector
        atmosphere.temperature[k] = T(288.15) - T(6.5) * (k-1) * T(0.5)  # Simple lapse rate
    end
    
    # Absorber profiles
    for k in 1:n_layers
        # Water vapor (H2O) - first absorber
        atmosphere.absorber[k, 1] = T(0.01) * exp(-atmosphere.pressure[k] / T(500))  # Simplified
        
        # Carbon dioxide (CO2) - constant
        atmosphere.absorber[k, 2] = T(400e-6)  # 400 ppm
        
        # Ozone (O3) - stratospheric peak
        if atmosphere.pressure[k] < 100
            atmosphere.absorber[k, 3] = T(5e-6)  # Stratospheric ozone
        else
            atmosphere.absorber[k, 3] = T(0.5e-6)  # Tropospheric ozone
        end
        
        # Other absorbers (simplified)
        atmosphere.absorber[k, 4] = T(320e-9)   # N2O
        atmosphere.absorber[k, 5] = T(1.8e-7)   # CO
        atmosphere.absorber[k, 6] = T(1.8e-6)   # CH4
    end
    
    # Altitude calculation (hydrostatic)
    R_dry = T(287.05)  # J/(kg·K)
    g = T(9.80665)     # m/s²
    
    atmosphere.level_altitude[n_layers+1] = location.altitude / 1000  # Convert to km
    for k in n_layers:-1:1
        # Hydrostatic equation
        scale_height = R_dry * atmosphere.temperature[k] / g
        dz = scale_height * log(pressure_levels[k+1] / pressure_levels[k])
        atmosphere.level_altitude[k] = atmosphere.level_altitude[k+1] + dz / 1000  # km
    end
    
    # Layer-mean altitudes
    for k in 1:n_layers
        atmosphere.altitude[k] = (atmosphere.level_altitude[k] + atmosphere.level_altitude[k+1]) / 2
    end
    
    return atmosphere
end

"""
    create_surface_properties(state::AbstractStateVector{T},
                             location::ObservationLocation{T}) where T

Create CRTM surface properties from model state at observation location.

# Arguments
- `state::AbstractStateVector{T}`: Model state vector
- `location::ObservationLocation{T}`: Observation location

# Returns
- `CRTMSurface{T}`: CRTM surface properties structure
"""
function create_surface_properties(state::AbstractStateVector{T},
                                   location::ObservationLocation{T}) where T
    
    surface = CRTMSurface{T}()
    
    # Determine surface type based on location
    if location.latitude < 60 && location.latitude > -60
        if location.longitude >= -180 && location.longitude <= 180
            # Simplified land/ocean classification
            if abs(location.longitude) < 30  # Rough Atlantic/African sector
                surface.surface_type = LAND_SURFACE
                surface.land_coverage = T(0.8)
                surface.water_coverage = T(0.2)
            else
                surface.surface_type = OCEAN_SURFACE
                surface.land_coverage = T(0.1)
                surface.water_coverage = T(0.9)
            end
        end
    else
        # Polar regions - consider ice coverage
        surface.surface_type = OCEAN_SURFACE
        surface.ice_coverage = T(0.5)
        surface.water_coverage = T(0.5)
    end
    
    # Extract surface temperature from model state
    # This would use proper interpolation from state vector
    surface.water_temperature = T(285.0)  # Simplified
    surface.land_temperature = T(283.0)   # Simplified
    
    # Wind speed from model state  
    # This would interpolate 10-meter winds
    surface.wind_speed = T(7.0)     # Simplified
    surface.wind_direction = T(225.0)  # Simplified
    
    # Other surface properties (would be from auxiliary data)
    surface.salinity = T(35.0)                    # Typical ocean salinity
    surface.vegetation_fraction = T(0.6)          # If land surface
    surface.soil_moisture_content = T(0.25)       # If land surface
    surface.lai = T(2.5)                         # Leaf area index
    
    return surface
end

"""
    create_sensor_geometry(observation::RadianceObservation{T}) where T

Create CRTM sensor geometry from satellite radiance observation.

# Arguments
- `observation::RadianceObservation{T}`: Satellite radiance observation

# Returns
- `CRTMGeometry{T}`: CRTM sensor geometry structure
"""
function create_sensor_geometry(observation::RadianceObservation{T}) where T
    
    return CRTMGeometry{T}(
        observation.satellite_zenith_angle,        # sensor_zenith_angle
        observation.satellite_azimuth_angle,       # sensor_azimuth_angle
        observation.solar_zenith_angle,            # solar_zenith_angle
        observation.solar_azimuth_angle,           # solar_azimuth_angle
        observation.scan_angle,                    # sensor_scan_angle
        observation.scan_position,                 # sensor_scan_position
        1,                                         # ifov (field of view)
        observation.location.latitude,             # latitude
        observation.location.longitude,            # longitude  
        observation.location.altitude              # surface_altitude
    )
end

# Jacobian computation functions (simplified implementations)

function compute_temperature_jacobian(layer::Int, channel::ChannelInfo{T},
                                     atmosphere::CRTMAtmosphere{T},
                                     geometry::CRTMGeometry{T}) where T
    # Simplified temperature Jacobian
    # Real implementation would use CRTM-computed Jacobians
    
    # Weighting function depends on channel frequency and pressure
    pressure = atmosphere.pressure[layer]
    
    if channel.central_frequency > 50  # High-peaking MW channels
        weight = exp(-(pressure - T(10))^2 / (2 * T(100)^2))  # Peak around 10 hPa
    else  # Surface-influenced channels
        weight = exp(-(pressure - T(850))^2 / (2 * T(200)^2))  # Peak around 850 hPa
    end
    
    # Zenith angle dependency
    zenith_factor = 1 / cos(deg2rad(geometry.sensor_zenith_angle))
    
    return weight * zenith_factor * T(0.1)  # Simplified scaling
end

function compute_humidity_jacobian(layer::Int, channel::ChannelInfo{T},
                                  atmosphere::CRTMAtmosphere{T}, 
                                  geometry::CRTMGeometry{T}) where T
    # Simplified humidity Jacobian for water vapor channels
    
    pressure = atmosphere.pressure[layer]
    
    # Humidity channels typically peak in mid-troposphere
    weight = exp(-(pressure - T(400))^2 / (2 * T(200)^2))  # Peak around 400 hPa
    
    # Only significant for humidity-sensitive channels
    if channel.central_frequency < 25 || channel.central_frequency > 180  # MW humidity channels
        weight *= T(0.05)
    else
        weight *= T(0.001)  # Minimal sensitivity
    end
    
    return weight * T(1.0)
end

function compute_surface_jacobian(channel::ChannelInfo{T},
                                 surface::CRTMSurface{T},
                                 geometry::CRTMGeometry{T}) where T
    # Surface temperature Jacobian for window channels
    
    # Strong for low-frequency MW and window IR channels
    if channel.central_frequency < 50  # Low-frequency MW
        surface_weight = T(0.8)
    elseif 900 < channel.central_frequency < 1100  # IR window
        surface_weight = T(0.9)
    else
        surface_weight = T(0.1)  # Atmospheric channels
    end
    
    # Zenith angle dependency (surface emissivity)
    emissivity_factor = 1 - T(0.3) * cos(deg2rad(geometry.sensor_zenith_angle))^2
    
    return surface_weight * emissivity_factor
end

# Instrument configuration creators

"""
    create_amsua_config() -> InstrumentConfig{Float64}

Create configuration for AMSU-A instrument.
"""
function create_amsua_config()
    sensor_info = SensorInfo(
        "amsua_n19",           # sensor_id
        "NOAA-19",             # platform_name  
        "AMSU-A",              # instrument_name
        :microwave,            # sensor_type
        15,                    # n_channels
        collect(1:15),         # channel_list
        [true, true, true, true, true, true, true, true, 
         true, true, true, true, true, true, true]  # use_flag
    )
    
    # AMSU-A channel characteristics
    frequencies = [23.8, 31.4, 50.3, 52.8, 53.596, 54.40, 54.94, 55.50,
                   57.290344, 57.290344, 57.290344, 57.290344, 57.290344, 89.0, 150.0]
    
    channel_info = [ChannelInfo{Float64}(
        i,                     # channel_number
        frequencies[i],        # central_frequency
        0.3,                   # noise_equivalent_delta_t
        true,                  # use_flag
        1.0,                   # error_variance
        0.0                    # bias_correction
    ) for i in 1:15]
    
    return InstrumentConfig{Float64}(
        sensor_info,
        channel_info,
        "amsua_n19.SpcCoeff.bin",    # coefficient_file
        nothing,                      # solar_irradiance_file
        :NoAerosolScattering,        # aerosol_option
        :NoCloudScattering          # cloud_option
    )
end

"""
    create_iasi_config() -> InstrumentConfig{Float64}

Create configuration for IASI instrument.
"""
function create_iasi_config()
    # IASI has 8461 channels, but typically only a subset is used
    selected_channels = [16, 29, 32, 35, 38, 41, 44, 47, 50, 55, 57, 59, 61, 63, 66, 70, 72, 74, 76, 78,
                        81, 83, 85, 87, 89, 92, 95, 97, 99, 101, 104, 106, 109, 111, 113, 116, 119, 122, 125, 128]
    
    sensor_info = SensorInfo(
        "iasi_metop-a",        # sensor_id
        "Metop-A",             # platform_name
        "IASI",                # instrument_name
        :infrared,             # sensor_type
        length(selected_channels),  # n_channels
        selected_channels,     # channel_list
        fill(true, length(selected_channels))  # use_flag
    )
    
    # IASI channels span 645-2760 cm⁻¹
    channel_info = [ChannelInfo{Float64}(
        ch,                    # channel_number
        645.0 + (ch-1) * 0.25, # central_frequency (cm⁻¹)
        0.5,                   # noise_equivalent_delta_t
        true,                  # use_flag
        2.0,                   # error_variance  
        0.0                    # bias_correction
    ) for ch in selected_channels]
    
    return InstrumentConfig{Float64}(
        sensor_info,
        channel_info,
        "iasi_metop-a.SpcCoeff.bin",     # coefficient_file
        nothing,                          # solar_irradiance_file
        :AerosolScattering,              # aerosol_option
        :CloudScattering                 # cloud_option
    )
end

end # module CRTMInterface